Stick Insect Climate PCA

Sorry for all of the packages. They grew and grew and grew, and I don’t want to refactor everything to reduce the dependencies.

source(file.path(r_path, "move_to_species.R"))

Error in file.path(r_path, "move_to_species.R") : 
  object 'r_path' not found

Read in the spreadsheet and take a look at the data.

# data path
loc_path <- here::here("data", "all species New_6-14-19.xlsx")

### read in spreadsheet
loc <- read_xlsx(loc_path) %>%
  janitor::clean_names() %>% 
  mutate(reproductive_mode = as.factor(reproductive_mode)) 

#get the number of individuals, and the sexuality counts per species
count_repro_mode <- loc %>% 
  group_by(genus, species, reproductive_mode) %>% 
  dplyr::count() %>% 
  mutate(genus_species = str_c(genus, species, sep = "_"),
         genus_species = str_replace_all(genus_species, " ", "_"),
         genus_species = str_replace_all(genus_species, "\\.", "")) %>% 
  ungroup() %>% 
  mutate(genus_species = fct_reorder(genus_species, n, sum)) %>% 
  ggplot(aes(x = genus_species, y = n, fill = reproductive_mode)) +
  geom_col() +
  coord_flip() + 
  theme_minimal()

count_repro_mode

Map

Plot a leaflet map of the localities. The leaflet map is interactive. You can click on the localities and a flag with some metadata will pop up!

#make locality shape file and assign WGS coord system
coord_points <- st_as_sf(loc, coords = c("longitude", "latitude"), 
                         crs = 4326, agr = "constant")

#use sourced plot_locs_leaflet script to plot localities
all_plot <- plot_locs_leaflet(loc, "reproductive_mode")

all_plot

# save the map
interactive_path <- here::here("output", "interactive_plots")
mapview::mapshot(
  all_plot,
  url = file.path(interactive_path, "all_species_map.html"),
  file = file.path(interactive_path, "all_species_map.pdf")
)

PCA-Genera

Climate Data

Read in the bioclim layers for analysis. I’m using all 19 for this preliminary exploration. I’m using CHELSA v1.2 data downloaded from their website. Plotting the first temperature layer to take a look at the data.

clim_files <- here::here("data", "climate") %>% 
  list.files(pattern = ".tif", full.names = TRUE)

w <- stack(clim_files)

#shorten the names of the bioclims
names(w) <- paste0(
  stringr::str_extract(names(w), "bio"), 
  stringr::str_extract(names(w), "\\d+$")
  )


ggplot() +
  ggspatial::layer_spatial(data = w[[1]]) +
  scale_fill_viridis_c(na.value = "transparent") +
  labs(fill = "Annual Avg Temp (C*10)") +
  theme_minimal()
  

PCA by locality

This is a PCA of the climate data extracted for each locality, rather than a PCA of the total climate space. This gives us a general idea of differences in environmental niche.

Run the pca and check out variable loadings and proportion of variance explained by components.

#extract data from chelsa for each locality. Making this into a data frame with columns labeled so the row labeling lines up after I remove the NAs.
#extract data from worldclim for each locality.
coords <- data.frame(latitude = loc$longitude, longitude = loc$latitude)

loc_clim <- dplyr::bind_cols(loc, raster::extract(w, coords, method = "simple", df = TRUE)) %>% 
  drop_na(bio1) %>% 
  dplyr::select(-ID)

#make a matrix of only bioclim values
clim_mat <- loc_clim[,grep("bio", names(loc_clim))] %>% as.matrix()

#run pca on climate variables
clim_pca <- prcomp(clim_mat, scale = TRUE)
summary_pca <- summary(clim_pca) #check out the components

#plot tables
summary_pca
knitr::kable(round(clim_pca$rotation[,1:3],3)) #Table of loading scores for the first 3 PCs.

Two plots: One plot of the PCA colored according to genus, with convex hulls surrounding the genera. It looks like this reflects a latitudinal gradient in temperature! You can interact with the PCA plot by clicking on points to view associated metadata. You can isolate the genus you want to view by double clicking the genus in the legend! You can also remove a genus by clicking on it once. There’s some other functionality you can explore in the toolbar at the top of the plot. The second plot is a PCA colored according to reproductive mode. It looks like asexual populations occupy slightly larger niche space, but both reproductive modes have a similar niche center.

#add pca results to loc_clim data frame
loc_clim <- data.frame(loc_clim, clim_pca$x)

#use sourced plot_clim_pca function to plot the pca results. args are the data set with species names and PC axis values and the pca summary
all_pca <- plot_clim_pca(loc_clim, summary_pca, factor = "genus")
all_pca

#use sourced plot_clim_pca function to plot the pca results. args are the data set with species names and PC axis values and the pca summary
repro_pca <-
  plot_clim_pca(loc_clim, summary_pca, factor = "reproductive_mode")
repro_pca


# save the plot colored by genus
htmlwidgets::saveWidget(
  all_pca,
  file.path(interactive_path, "all_species_pca_genus.html"),
  selfcontained = TRUE
)

# save the plot colored by reproductive mode
htmlwidgets::saveWidget(
  repro_pca,
  file.path(interactive_path, "all_species_pca_repro.html"),
  selfcontained = TRUE
)

Examining whether asexual populations occupy more unstable climates than sexual populations. Only using species with multiple sexual and asexual populations. Asexual pops tend to occupy more stable temperature environments, but less stable precipitation environments. We’re estimating stability using the method presented by Owens and Guralnick 2019- climateStability: AN R PACKAGE TO ESTIMATE CLIMATE STABILITY FROM TIME-SLICE CLIMATOLOGIES.

### Creating temperature and precipitation stability layers
### Using bio1 (average temp) and bio12 (average precip)
### 2.5 arcmin resolution, already cropped to NZ to speed up computation time


#time slices are calculated as the difference between the midpoints of the two time periods the climate layers were calculated for (e.g. midpoint of LH = (4.2 ka - 0.3 ka) / 2 + 0.3 ka = 2.25, midpoint of MH = (8.326 ka - 4.2 ka) / 2 + 4.2 = 6.263. time_slice = 6.263 - 2.25 = 4.013)
time_slices <- c(2550, 4013, 6037, 1500, 2050, 5150)

precip_path <-
  here::here("data",
       "climate",
       "paleoclim_late_pleistocene",
       "precipitation")

precip_deviation <-
  climateStability::deviationThroughTime(variableDirectory = precip_path, timeSlicePeriod = time_slices)

precip_stability <- (1 / precip_deviation)
precip_stability_scaled <- climateStability::rescale0to1(precip_stability)

temp_path <-
  here::here("data",
       "climate",
       "paleoclim_late_pleistocene",
       "temperature")

temp_deviation <-
  climateStability::deviationThroughTime(variableDirectory = temp_path, timeSlicePeriod = time_slices)

temp_stability <- (1 / temp_deviation)
temp_stability_scaled <- climateStability::rescale0to1(temp_stability)

overall_stability_scaled <- precip_stability_scaled * temp_stability_scaled

temp_stability_map <- ggplot() +
  ggspatial::layer_spatial(data = temp_stability_scaled) +
  labs(title = "Temperature stability") +
  scale_fill_viridis_c(na.value = "transparent") +
  theme_minimal()

precip_stability_map <- ggplot() +
  ggspatial::layer_spatial(data = precip_stability_scaled) +
  labs(title = "Precipitation stability") +
  scale_fill_viridis_c(na.value = "transparent") +
  theme_minimal()

overall_stability_map <- ggplot() +
  ggspatial::layer_spatial(data = overall_stability_scaled) +
  labs(title = "Overall stability") +
  scale_fill_viridis_c(na.value = "transparent") +
  theme_minimal()

temp_stability_map
precip_stability_map
overall_stability_map

# write maps to file
map_path <- here::here("output", "maps")

ggsave("temp_stability.png",
       plot = temp_stability_map,
       device = "png",
       path = map_path,
       dpi = "retina")

ggsave("precip_stability.png",
       plot = precip_stability_map,
       device = "png",
       path = map_path,
       dpi = "retina")

ggsave("overall_stability.png",
       plot = overall_stability_map,
       device = "png",
       path = map_path,
       dpi = "retina")

Plot stability across species.

# filter for relevant species and asexual reproductive mode
asexual_locs <- loc %>%
  mutate(lat_long = str_c(latitude, longitude)) %>%
  filter(
    reproductive_mode == "asexual",
    species == "horridus" |
      species == "jucundum" |
      species == "hookeri" |
      species == "annulata" |
      species == "ovobessus" |
      species == "huttoni",
    !duplicated(lat_long)
  ) %>%
  dplyr::select(longitude, latitude)

# filter for relevant species and sexual reproductive mode
sexual_locs <- loc %>%
  mutate(lat_long = str_c(latitude, longitude)) %>%
  filter(
    reproductive_mode == "sexual",
    species == "horridus" |
      species == "jucundum" |
      species == "hookeri" |
      species == "annulata" |
      species == "ovobessus" |
      species == "huttoni",
    !duplicated(lat_long)
  ) %>%
  dplyr::select(longitude, latitude)

# extract preciptitation values and bind into a new dataframe
precip_asexual <-
  raster::extract(precip_stability_scaled, asexual_locs) %>%
  enframe(name = NULL, value = "precip_stability_scaled") %>%
  mutate(reproductive_mode = "asexual")

precip_sexual <-
  raster::extract(precip_stability_scaled, sexual_locs) %>%
  enframe(name = NULL, value = "precip_stability_scaled") %>%
  mutate(reproductive_mode = "sexual")

precip_df <- bind_rows(precip_asexual, precip_sexual)

# plot precipitation stability
precip_stability_plot <- ggplot(
  data = precip_df,
  aes(x = reproductive_mode, y = precip_stability_scaled, fill = reproductive_mode)
) +
  geom_violin(width = 0.8) +
  geom_boxplot(width = 0.1,
               color = "gray",
               fill = "transparent") +
  scale_fill_viridis_d(option = "magma") +
  theme_dark()

precip_stability_plot

# extract temperature values and bind into a new data frame
temp_asexual <-
  raster::extract(temp_stability_scaled, asexual_locs) %>%
  enframe(name = NULL, value = "temp_stability_scaled") %>%
  mutate(reproductive_mode = "asexual")

temp_sexual <-
  raster::extract(temp_stability_scaled, sexual_locs) %>%
  enframe(name = NULL, value = "temp_stability_scaled") %>%
  mutate(reproductive_mode = "sexual")

temp_df <- bind_rows(temp_asexual, temp_sexual)

# plot temperature stability
temp_stability_plot <- ggplot(data = temp_df,
       aes(x = reproductive_mode, y = temp_stability_scaled, fill = reproductive_mode)) +
  geom_violin(width = 0.8) +
  geom_boxplot(width = 0.1,
               color = "gray",
               fill = "transparent") +
  scale_fill_viridis_d(option = "magma") +
  theme_dark()

temp_stability_plot

# extract overall stability values and bind into a dataframe
overall_asexual <-
  raster::extract(overall_stability_scaled, asexual_locs) %>%
  enframe(name = NULL, value = "overall_stability_scaled") %>%
  mutate(reproductive_mode = "asexual")

overall_sexual <-
  raster::extract(overall_stability_scaled, sexual_locs) %>%
  enframe(name = NULL, value = "overall_stability_scaled") %>%
  mutate(reproductive_mode = "sexual")

overall_df <- bind_rows(overall_asexual, overall_sexual)

# plot overall stability
overall_stability_plot <- ggplot(
  data = overall_df,
  aes(x = reproductive_mode, y = overall_stability_scaled, fill = reproductive_mode)
) +
  geom_violin(width = 0.8) +
  geom_boxplot(width = 0.1,
               color = "gray",
               fill = "transparent") +
  scale_fill_viridis_d(option = "magma") +
  theme_dark()

overall_stability_plot

# save plots
plot_path <- here::here("output", "plots")
ggsave("temp_stability_plot.png",
       plot = temp_stability_plot,
       device = "png",
       path = plot_path,
       dpi = "retina")

ggsave("precip_stability_plot.png",
       plot = precip_stability_plot,
       device = "png",
       path = plot_path,
       dpi = "retina")

ggsave("overall_stability_plot.png",
       plot = overall_stability_plot,
       device = "png",
       path = plot_path,
       dpi = "retina")

PCA-Species

These are PCAs of environmental space for species within genera. Each climate PCA is of localities for a single genus, colored by species. I’m doing this even for genera with one species, so it’s easy to see if certain localities group together.

Acanthoxyla

#source function to conduct a PCA on individual species
summary_list_acan <- species_pca_fun(loc_clim, "acanthoxyla")
#plot
acan_plot <- plot_clim_pca(summary_list_acan$loc_clim, summary_list_acan$summary_pca, "reproductive_mode")

acan_plot

# save pca plot
htmlwidgets::saveWidget(acan_plot, file.path(interactive_path, "acanthoxyla_pca.html"), selfcontained = TRUE)

#filter localities for the focal genus
acan_loc <- loc %>% 
  filter(genus == "acanthoxyla")
  
#use sourced plot_locs_leaflet script to plot localities
acan_map <- plot_locs_leaflet(acan_loc, "reproductive_mode")

acan_map

# save the map
mapview::mapshot(acan_map, url = file.path(interactive_path, "acan_map.html"), file = file.path(interactive_path,"acan_map.pdf"))

summary_list_acan$summary_pca
loadings_acan <- summary_list_acan$summary_pca$rotation
knitr::kable(round(loadings_acan[,1:3],3)) #Table of loading scores for the first 3 PCs. 

Argosarchus

# conduct pca
summary_list_argo <- species_pca_fun(loc_clim, "argosarchus")

# plot
argo_plot <-
  plot_clim_pca(summary_list_argo$loc_clim,
                summary_list_argo$summary_pca,
                factor = "reproductive_mode")

Ignoring unknown aesthetics: text

argo_plot

htmlwidgets::saveWidget(argo_plot,
                        file.path(interactive_path, "argosarchus_pca.html"),
                        selfcontained = TRUE)

#filter localities for the focal genus
argo_loc <- loc %>%
  filter(genus == "argosarchus")

#use sourced plot_locs_leaflet script to plot localities
argo_map <- plot_locs_leaflet(argo_loc, "reproductive_mode")

Assuming "longitude" and "latitude" are longitude and latitude, respectively

argo_map

mapview::mapshot(
  argo_map,
  url = file.path(interactive_path, "argo_map.html"),
  file = file.path(interactive_path, "argo_map.pdf")
)

summary_list_argo$summary_pca

Importance of components:
                          PC1    PC2    PC3    PC4     PC5     PC6
Standard deviation     2.8521 2.4343 1.6374 1.0915 0.76841 0.51224
Proportion of Variance 0.4281 0.3119 0.1411 0.0627 0.03108 0.01381
Cumulative Proportion  0.4281 0.7400 0.8811 0.9438 0.97491 0.98872
                           PC7    PC8    PC9    PC10    PC11    PC12
Standard deviation     0.35867 0.2223 0.1072 0.09909 0.08123 0.05326
Proportion of Variance 0.00677 0.0026 0.0006 0.00052 0.00035 0.00015
Cumulative Proportion  0.99549 0.9981 0.9987 0.99921 0.99956 0.99971
                          PC13    PC14    PC15    PC16    PC17
Standard deviation     0.04063 0.03890 0.03293 0.02443 0.01689
Proportion of Variance 0.00009 0.00008 0.00006 0.00003 0.00002
Cumulative Proportion  0.99980 0.99988 0.99993 0.99996 0.99998
                          PC18    PC19
Standard deviation     0.01643 0.01112
Proportion of Variance 0.00001 0.00001
Cumulative Proportion  0.99999 1.00000

loadings_argo <- summary_list_argo$summary_pca$rotation
knitr::kable(round(loadings_argo[, 1:3], 3)) #Table of loading scores for the first 3 PCs. 

	PC1	PC2	PC3
bio1	-0.266	-0.227	-0.182
bio10	-0.242	-0.207	-0.292
bio11	-0.280	-0.228	-0.110
bio12	-0.225	0.313	-0.036
bio13	-0.257	0.273	-0.023
bio14	-0.202	0.330	-0.026
bio15	-0.141	-0.175	0.020
bio16	-0.254	0.275	-0.033
bio17	-0.197	0.335	-0.030
bio18	-0.175	0.347	0.014
bio19	-0.276	0.219	-0.090
bio2	0.266	0.149	-0.319
bio3	0.268	0.109	-0.277
bio4	0.227	0.162	-0.310
bio5	-0.113	-0.143	-0.516
bio6	-0.292	-0.224	-0.008
bio7	0.258	0.163	-0.324
bio8	-0.098	-0.073	0.249
bio9	-0.191	-0.155	-0.385

Now I’m going to to environmental niche factor analysis between sexual and asexual populations within the species.

# create path for spreadsheets
spread_path <- here::here("output", "spreadsheets")

#get background env't for the species
ahor_bg_env <- bg_env_crop(argo_loc, 
                           species = "horridus",
                           environment = w, 
                           buffer = 0.5)

#enfa for the sexual species
ahor_sexual_enfa <- enfa_calc_fun(locs = argo_loc, 
                                  species = "horridus", 
                                  reproductive_mode = "sexual", 
                                  mask_raster = ahor_bg_env)

#enfa for the asexual species
ahor_asexual_enfa <- enfa_calc_fun(locs = argo_loc, 
                                   species = "horridus", 
                                   reproductive_mode = "asexual", 
                                   mask_raster = ahor_bg_env)


#plot the marginality scores
ahor_marginality <- marginality_lollipop(sex_marg = ahor_sexual_enfa$m, 
                    asex_marg = ahor_asexual_enfa$m,
                    full_species_name = "Argosarchus horridus")

# write scores to csvs
readr::write_csv(ahor_asexual_enfa$m %>% enframe(), 
                 file.path(spread_path, "ahor_asexual_marginality_score.csv"))

readr::write_csv(ahor_sexual_enfa$m %>% enframe(), 
                 file.path(spread_path, "ahor_sexual_marginality_score.csv"))

ggsave("ahor_marginality_lollipop.png",
       plot = ahor_marginality,
       device = "png",
       path = plot_path,
       dpi = "retina")

A couple different ways to visualize the environmental variation. 1) Scatterplot visualizations of marginality vs axis 1 of the specialization with the labels removed (they make things indiscernable). Red = occupied e-space. Gray = Background e-space. The yellow dot indicates the mean marginality (it’s not the value that is on the lollipop plot, so don’t let that confuse you). 2) ENFA histogram visualizations of marginality and specialization axes. 3) PCA of total e-space with colors corresponding to sexual vs. asexual populations.

### 1) ENFA scatterplot
#access the relevant values for plotting
ahor_asexual_df <- ahor_asexual_enfa$li %>% 
  as_tibble() %>% 
  bind_cols(pr = ahor_asexual_enfa$pr)

readr::write_csv(ahor_asexual_df, 
                 file.path(spread_path, "ahor_asexual_marginality_df.csv"))

ahor_sexual_df <- ahor_sexual_enfa$li %>% 
  as_tibble() %>% 
  bind_cols(pr = ahor_sexual_enfa$pr)

readr::write_csv(ahor_sexual_df, 
                 file.path(spread_path, "ahor_sexual_marginality_df.csv"))

#asexual
ahor_enfa_spec_asexual <- enfa_hex_plot(ahor_asexual_df, marg = Mar, spec = Spe1, repro_mode = "Asexual")

ahor_enfa_spec_asexual

ggsave("ahor_enfa_spec_asexual.png",
       plot = ahor_enfa_spec_asexual,
       device = "png",
       path = plot_path,
       dpi = "retina")

#sexual
ahor_enfa_spec_sexual <- enfa_hex_plot(ahor_sexual_df, marg = Mar, spec = Spe1, repro_mode = "Sexual")

ahor_enfa_spec_sexual

ggsave("ahor_enfa_spec_sexual.png",
       plot = ahor_enfa_spec_sexual,
       device = "png",
       path = plot_path,
       dpi = "retina")

### 2) ENFA histogram
## asexual
hist(ahor_asexual_enfa)
title(main = "Asexual", adj = 0.7, line = -12)

# write to file
png(filename = file.path(plot_path, "ahor_asexual_enfa_hist.png"))
hist(ahor_asexual_enfa)

title(main = "Asexual", adj = 0.7, line = -12)
dev.off()

quartz_off_screen 
                2 

## sexual
hist(ahor_sexual_enfa)
title(main = "Sexual", adj = 0.7, line = -12)

# write to file
png(filename = file.path(plot_path, "ahor_sexual_enfa_hist.png"))
hist(ahor_sexual_enfa)

title(main = "Sexual", adj = 0.7, line = -12)
dev.off()

quartz_off_screen 
                2 

### 3) PCA of total e-space
ahor_total_pca <- total_climate_pca_plot(bg_env = ahor_bg_env, locs = argo_loc, genus = "Argosarchus", species = "horridus")

ahor_total_pca

ggsave("ahor_total_pca.png",
       plot = ahor_total_pca,
       device = "png",
       path = plot_path,
       dpi = "retina")

Trying out a repeat of the analyses with reduced environmental space. Prioritizing variables that will limit their distribution (i.e. variables that represent the extremes). After correlation analysis, we’re left with BIO6, BIO13, BIO14, BIO16

# PCA of background e-space. Resulting list is a correlation heatmap (cor_heatmap), a tibble of the rasters with correlation > the cutoff (default is 0.75), and a tibble of all pairwise correlations
ahor_pca <- raster_correlation(raster_stack = ahor_bg_env)

Warning in for (b in 1:raster@data@nlayers) { :
  closing unused connection 5 (/private/var/folders/vf/3z__6c_11bb7b95dtcms1m0c0000gp/T/Rtmp6hXykz/raster/r_tmp_2020-03-06_091749_61578_87059.gri)

ahor_pca$cor_heatmap
ahor_pca$top_cor %>% knitr::kable()

Var1	Var2	value
bio5	bio9	0.8044233
bio11	bio9	0.7868951
bio11	bio5	0.8900737
bio6	bio9	0.7577603
bio6	bio5	0.8347273
bio6	bio11	0.9932021
bio1	bio9	0.7957997
bio1	bio5	0.9366271
bio1	bio11	0.9917017
bio1	bio6	0.9736459
bio10	bio9	0.8125594
bio10	bio5	0.9720882
bio10	bio11	0.9672682
bio10	bio6	0.9367587
bio10	bio1	0.9910415
bio12	bio19	0.9269863
bio13	bio19	0.9488991
bio13	bio12	0.9924923
bio16	bio19	0.9527998
bio16	bio12	0.9922873
bio16	bio13	0.9990530
bio18	bio19	0.8453521
bio18	bio12	0.9789417
bio18	bio13	0.9609296
bio18	bio16	0.9575710
bio14	bio19	0.8942461
bio14	bio12	0.9909216
bio14	bio13	0.9755415
bio14	bio16	0.9742782
bio14	bio18	0.9868791
bio17	bio19	0.8939912
bio17	bio12	0.9921150
bio17	bio13	0.9758882
bio17	bio16	0.9755773
bio17	bio18	0.9838595
bio17	bio14	0.9974118
bio2	bio3	0.9186390
bio2	bio4	0.8986994
bio7	bio3	0.8312469
bio7	bio4	0.9628226
bio7	bio2	0.9802946

ggsave("ahor_pca_corr.png",
       plot = ahor_pca$cor_heatmap,
       device = "png",
       path = plot_path,
       dpi = "retina")

Saving 5.97 x 3.69 in image

# create path for spreadsheets
spread_path <- here::here("output", "spreadsheets")

# write correlation data frame to file
readr::write_csv(ahor_pca$top_cor, file.path(spread_path, "ahor_top_cor.csv"))

Repeat the analysis with the reduced data set. The background environment is 0.5 degrees, a ballpark dispersal limitation for all stick insect species in this study.

w_ahor <- raster::subset(w, c("bio6", "bio13", "bio14", "bio16"))

#get background env't for the species
ahor_bg_env_subset <- bg_env_crop(argo_loc, 
                           species = "horridus",
                           environment = w_ahor, 
                           buffer = 0.5)

#enfa for the sexual species
ahor_sexual_enfa_subset <- enfa_calc_fun(locs = argo_loc, 
                                  species = "horridus", 
                                  reproductive_mode = "sexual", 
                                  mask_raster = ahor_bg_env_subset)

#enfa for the asexual species
ahor_asexual_enfa_subset <- enfa_calc_fun(locs = argo_loc, 
                                   species = "horridus", 
                                   reproductive_mode = "asexual", 
                                   mask_raster = ahor_bg_env_subset)

# write marginality scores to csv
readr::write_csv(
  ahor_asexual_enfa_subset$m %>% enframe(name = NULL, value = "marginality"),
  file.path(spread_path, "ahor_subset_asexual_marginality_score.csv")
)

readr::write_csv(
  ahor_sexual_enfa_subset$m %>% enframe(name = NULL, value = "marginality"),
  file.path(spread_path, "ahor_subset_sexual_marginality_score.csv")
)

#plot the marginality scores
ahor_subset_marginality <-
  marginality_lollipop(
    sex_marg = ahor_sexual_enfa_subset$m,
    asex_marg = ahor_asexual_enfa_subset$m,
    full_species_name = "Argosarchus horridus"
  )

ahor_subset_marginality

ggsave("ahor_subset_marginality.png",
       plot = ahor_subset_marginality,
       device = "png",
       path = plot_path,
       dpi = "retina")

Saving 5.97 x 3.69 in image

Visualize with reduced data set

### 1) ENFA scatterplot
#access the relevant values for plotting
ahor_asexual_df_subset <- ahor_asexual_enfa_subset$li %>%
  as_tibble() %>%
  bind_cols(pr = ahor_asexual_enfa_subset$pr)

readr::write_csv(
  ahor_asexual_df_subset,
  file.path(spread_path, "ahor_subset_asexual_marginality_df.csv")
)

ahor_sexual_df_subset <- ahor_sexual_enfa_subset$li %>%
  as_tibble() %>%
  bind_cols(pr = ahor_sexual_enfa_subset$pr)

readr::write_csv(
  ahor_sexual_df_subset,
  file.path(spread_path, "ahor_subset_sexual_marginality_df.csv")
)

#asexual. Jesus these variable names are getting long
ahor_subset_enfa_spec_asexual <-
  enfa_hex_plot(
    ahor_asexual_df_subset,
    marg = Mar,
    spec = Spe1,
    repro_mode = "Asexual"
  )

ahor_subset_enfa_spec_asexual

ggsave(
  "ahor_subset_enfa_spec_asexual.png",
  plot = ahor_subset_enfa_spec_asexual,
  device = "png",
  path = plot_path,
  dpi = "retina"
)

#sexual
ahor_subset_enfa_spec_sexual <-
  enfa_hex_plot(
    ahor_sexual_df_subset,
    marg = Mar,
    spec = Spe1,
    repro_mode = "Sexual"
  )

ahor_subset_enfa_spec_sexual

ggsave(
  "ahor_subset_enfa_spec_sexual.png",
  plot = ahor_subset_enfa_spec_sexual,
  device = "png",
  path = plot_path,
  dpi = "retina"
)

### 2) ENFA histogram
# asexual
hist(ahor_asexual_enfa_subset)
title(main = "Asexual", adj = 0.7, line = -12)

png(filename = file.path(plot_path, "ahor_subset_asexual_enfa_hist.png"))
hist(ahor_asexual_enfa_subset)

title(main = "Asexual", adj = 0.7, line = -12)
dev.off()

quartz_off_screen 
                2 

# sexual
hist(ahor_sexual_enfa_subset)
title(main = "Sexual", adj = 0.7, line = -12)

png(filename = file.path(plot_path, "ahor_subset_sexual_enfa_hist.png"))
hist(ahor_sexual_enfa_subset)

title(main = "Sexual", adj = 0.7, line = -12)
dev.off()

quartz_off_screen 
                2 

### 3) PCA of total e-space
ahor_subset_total_pca <-
  total_climate_pca_plot(
    bg_env = ahor_bg_env_subset,
    locs = argo_loc,
    genus = "Argosarchus",
    species = "horridus"
  )

ahor_subset_total_pca

ggsave(
  "ahor_subset_total_pca.png",
  plot = ahor_subset_total_pca,
  device = "png",
  path = plot_path,
  dpi = "retina"
)

# output object path
obj_path <- here::here("output", "objects")

# save enfa objects and remove them from the environment. They take up a lot of memory.
saveRDS(ahor_sexual_enfa, file = file.path(obj_path, "ahor_sexual_enfa.RDS"))
rm(ahor_sexual_enfa)

saveRDS(ahor_asexual_enfa, file = file.path(obj_path, "ahor_asexual_enfa.RDS"))
rm(ahor_asexual_enfa)

saveRDS(ahor_sexual_enfa_subset,
        file = file.path(obj_path, "ahor_subset_sexual_enfa.RDS"))
rm(ahor_sexual_enfa_subset)

saveRDS(ahor_asexual_enfa_subset,
        file = file.path(obj_path, "ahor_subset_asexual_enfa.RDS"))
rm(ahor_asexual_enfa_subset)

We’re also interested in seeing if asexual populations live in more unstable climatic areas relative to sexual populations.

argo_locs_asexual <- argo_loc %>% 
  mutate(lat_long = str_c(latitude, longitude, sep = "_")) %>% 
  filter(reproductive_mode == "asexual",
         !duplicated(lat_long)) %>% 
  dplyr::select(longitude, latitude)

argo_locs_sexual <- argo_loc %>% 
  mutate(lat_long = str_c(latitude, longitude, sep = "_")) %>% 
  filter(reproductive_mode == "sexual",
         !duplicated(lat_long)) %>% 
  dplyr::select(longitude, latitude)

precip_asexual_ahor <- raster::extract(precip_stability_scaled, argo_locs_asexual) %>% 
  enframe(name = NULL, value = "precip_stability_scaled") %>% 
  mutate(reproductive_mode = "asexual")

precip_sexual_ahor <- raster::extract(precip_stability_scaled, argo_locs_sexual) %>% 
  enframe(name = NULL, value = "precip_stability_scaled") %>% 
  mutate(reproductive_mode = "sexual")

precip_df_ahor <- bind_rows(precip_asexual_ahor, precip_sexual_ahor)

readr::write_csv(precip_df_ahor, 
                 file.path(spread_path, "ahor_precip_stability_df.csv"))

ahor_precip_stability_plot <- ggplot(data = precip_df_ahor, aes(x = reproductive_mode, y = precip_stability_scaled, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

ahor_precip_stability_plot

ggsave(
  "ahor_precip_stability.png",
  plot = ahor_precip_stability_plot,
  device = "png",
  path = plot_path,
  dpi = "retina"
)

temp_asexual_ahor <- raster::extract(temp_stability_scaled, argo_locs_asexual) %>% 
  enframe(name = NULL, value = "temp_stability_scaled") %>% 
  mutate(reproductive_mode = "asexual")

temp_sexual_ahor <- raster::extract(temp_stability_scaled, argo_locs_sexual) %>% 
  enframe(name = NULL, value = "temp_stability_scaled") %>% 
  mutate(reproductive_mode = "sexual")

temp_df_ahor <- bind_rows(temp_asexual_ahor, temp_sexual_ahor)

readr::write_csv(temp_df_ahor, 
                 file.path(spread_path, "ahor_temp_stability_df.csv"))

ahor_temp_stability_plot <- ggplot(data = temp_df_ahor, aes(x = reproductive_mode, y = temp_stability_scaled, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

ahor_temp_stability_plot

ggsave(
  "ahor_temp_stability.png",
  plot = ahor_temp_stability_plot,
  device = "png",
  path = plot_path,
  dpi = "retina"
)

overall_asexual_ahor <- raster::extract(overall_stability_scaled, argo_locs_asexual) %>% 
  enframe(name = NULL, value = "overall_stability_scaled") %>% 
  mutate(reproductive_mode = "asexual")

overall_sexual_ahor <- raster::extract(overall_stability_scaled, argo_locs_sexual) %>% 
  enframe(name = NULL, value = "overall_stability_scaled") %>% 
  mutate(reproductive_mode = "sexual")

overall_df_ahor <- bind_rows(overall_asexual_ahor, overall_sexual_ahor)

readr::write_csv(overall_df_ahor, 
                 file.path(spread_path, "ahor_overall_stability_df.csv"))


ahor_overall_stability_plot <- ggplot(data = overall_df_ahor, aes(x = reproductive_mode, y = overall_stability_scaled, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

ahor_overall_stability_plot

ggsave(
  "ahor_overall_stability.png",
  plot = ahor_overall_stability_plot,
  device = "png",
  path = plot_path,
  dpi = "retina"
)

Put all output into species-specific subfolders.

# move interactive
move_to_species(in_path = interactive_path,
                out_path = ahor_out_int_path,
                pattern = "ahor")

Asteliaphasma

#pca
summary_list_aste <- species_pca_fun(loc_clim, "asteliaphasma")
#plot
aste_plot <- plot_clim_pca(summary_list_aste$loc_clim, summary_list_aste$summary_pca, factor = "reproductive_mode")
aste_plot

#if selfcontained = TRUE, you can remove the folder that gets added alongside the plot. 
htmlwidgets::saveWidget(aste_plot, paste0(getwd(), "/plots/repro_mode_plots/asteliaphasma_pca.html"), selfcontained = TRUE)

#filter localities for the focal genus
aste_loc <- loc %>% 
  filter(genus == "asteliaphasma")
  
#use sourced plot_locs_leaflet script to plot localities
aste_map <- plot_locs_leaflet(aste_loc, "reproductive_mode")

aste_map

#in case I want to save the map somewhere
#mapview::mapshot(aste_map, url = paste0(getwd(), "/plots/repro_mode_plots/aste_map.html"), file = paste0(getwd(), "/plots/repro_mode_plots/aste_map.pdf"))

summary_list_aste$summary_pca
loadings_aste <- summary_list_aste$summary_pca$rotation
knitr::kable(round(loadings_aste[,1:3],3)) #Table of loading scores for the first 3 PCs. 

Now I’m going to to environmental niche factor analysis between sexual and asexual populations within the species.

#get background env't for the species
ajuc_bg_env <- bg_env_crop(aste_loc, 
                           species = "jucundum",
                           environment = w, 
                           buffer = 0.5)

#enfa for the sexual species
ajuc_sexual_enfa <- enfa_calc_fun(locs = aste_loc, 
                                  species = "jucundum", 
                                  reproductive_mode = "sexual", 
                                  mask_raster = ajuc_bg_env)

#enfa for the asexual species
ajuc_asexual_enfa <- enfa_calc_fun(locs = aste_loc, 
                                   species = "jucundum", 
                                   reproductive_mode = "asexual", 
                                   mask_raster = ajuc_bg_env)


#plot the marginality scores
marginality_lollipop(sex_marg = ajuc_sexual_enfa$m, 
                    asex_marg = ajuc_asexual_enfa$m,
                    full_species_name = "Asteliaphasma jucundum")

A couple different ways to visualize the environmental variation. 1) Scatterplot visualizations of marginality vs axis 1 of the specialization with the labels removed (they make things indiscernable). Red = occupied e-space. Gray = Background e-space. 2) ENFA histogram visualizations of marginality and specialization axes. 3) PCA of total e-space with colors corresponding to sexual vs. asexual populations.

### 1) ENFA scatterplot
#access the relevant values for plotting
ajuc_asexual_df <- ajuc_asexual_enfa$li %>% 
  as_tibble() %>% 
  bind_cols(pr = ajuc_asexual_enfa$pr)


ajuc_sexual_df <- ajuc_sexual_enfa$li %>% 
  as_tibble() %>% 
  bind_cols(pr = ajuc_sexual_enfa$pr)


#asexual
enfa_hex_plot(ajuc_asexual_df, marg = Mar, spec = Spe1, repro_mode = "Asexual")


#sexual
enfa_hex_plot(ajuc_sexual_df, marg = Mar, spec = Spe1, repro_mode = "Sexual")


### 2) ENFA histogram
#asexual
hist(ajuc_asexual_enfa)
title(main = "Asexual", adj = 0.7, line = -12)

#sexual
hist(ajuc_sexual_enfa)
title(main = "Sexual", adj = 0.7, line = -12)

### 3) PCA of total e-space
total_climate_pca_plot(bg_env = ajuc_bg_env, locs = aste_loc, genus = "Asteliophasma", species = "jucundum")

Trying out a repeat of the analyses with reduced environmental space. Prioritizing variables that will limit their distribution (i.e. variables that represent the extremes). After correlation analysis, we’re left with BIO5, BIO6, BIO14, BIO17.

#PCA of background e-space. Resulting list is a correlation heatmap (cor_heatmap), a tibble of the rasters with correlation > the cutoff (default is 0.75), and a tibble of all pairwise correlations
ajuc_pca <- raster_correlation(raster_stack = ajuc_bg_env)
ajuc_pca$cor_heatmap
ajuc_pca$top_cor %>% knitr::kable()

Repeat the analysis with the reduced data set.

w_ajuc <- raster::subset(w, c("bio5", "bio6", "bio14", "bio17"))

#get background env't for the species
ajuc_bg_env_subset <- bg_env_crop(aste_loc, 
                           species = "jucundum",
                           environment = w_ajuc, 
                           buffer = 0.5)

#enfa for the sexual species
ajuc_sexual_enfa_subset <- enfa_calc_fun(locs = aste_loc, 
                                  species = "jucundum", 
                                  reproductive_mode = "sexual", 
                                  mask_raster = ajuc_bg_env_subset)

#enfa for the asexual species
ajuc_asexual_enfa_subset <- enfa_calc_fun(locs = aste_loc, 
                                   species = "jucundum", 
                                   reproductive_mode = "asexual", 
                                   mask_raster = ajuc_bg_env_subset)


#plot the marginality scores
marginality_lollipop(sex_marg = ajuc_sexual_enfa_subset$m, 
                    asex_marg = ajuc_asexual_enfa_subset$m,
                    full_species_name = "Asteliaphasma jucundum")

Visualize with reduced data set

### 1) ENFA scatterplot
#access the relevant values for plotting
ajuc_asexual_df_subset <- ajuc_asexual_enfa_subset$li %>% 
  as_tibble() %>% 
  bind_cols(pr = ajuc_asexual_enfa_subset$pr)


ajuc_sexual_df_subset <- ajuc_sexual_enfa_subset$li %>% 
  as_tibble() %>% 
  bind_cols(pr = ajuc_sexual_enfa_subset$pr)


#asexual
enfa_hex_plot(ajuc_asexual_df_subset, marg = Mar, spec = Spe1, repro_mode = "Asexual")


#sexual
enfa_hex_plot(ajuc_sexual_df_subset, marg = Mar, spec = Spe1, repro_mode = "Sexual")


### 2) ENFA histogram
#asexual
hist(ajuc_asexual_enfa_subset)
title(main = "Asexual", adj = 0.7, line = -12)

#sexual
hist(ajuc_sexual_enfa_subset)
title(main = "Sexual", adj = 0.7, line = -12)

### 3) PCA of total e-space
total_climate_pca_plot(bg_env = ajuc_bg_env_subset, locs = aste_loc, genus = "Asteliaphasma", species = "jucundum")

We’re also interested in seeing if asexual populations live in more unstable climatic areas relative to sexual populations.

aste_locs_asexual <- aste_loc %>% 
  mutate(lat_long = str_c(latitude, longitude, sep = "_")) %>% 
  filter(reproductive_mode == "asexual",
         !duplicated(lat_long)) %>% 
  dplyr::select(longitude, latitude)

aste_locs_sexual <- aste_loc %>% 
  mutate(lat_long = str_c(latitude, longitude, sep = "_")) %>% 
  filter(reproductive_mode == "sexual",
         !duplicated(lat_long)) %>% 
  dplyr::select(longitude, latitude)

precip_asexual_ajuc <- raster::extract(precip_stability_scaled, aste_locs_asexual) %>% 
  enframe(name = NULL, value = "precip_stability_scaled") %>% 
  mutate(reproductive_mode = "asexual")

precip_sexual_ajuc <- raster::extract(precip_stability_scaled, aste_locs_sexual) %>% 
  enframe(name = NULL, value = "precip_stability_scaled") %>% 
  mutate(reproductive_mode = "sexual")

precip_df_ajuc <- bind_rows(precip_asexual_ajuc, precip_sexual_ajuc)

ggplot(data = precip_df_ajuc, aes(x = reproductive_mode, y = precip_stability_scaled, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()


temp_asexual_ajuc <- raster::extract(temp_stability_scaled, aste_locs_asexual) %>% 
  enframe(name = NULL, value = "temp_stability_scaled") %>% 
  mutate(reproductive_mode = "asexual")

temp_sexual_ajuc <- raster::extract(temp_stability_scaled, aste_locs_sexual) %>% 
  enframe(name = NULL, value = "temp_stability_scaled") %>% 
  mutate(reproductive_mode = "sexual")

temp_df_ajuc <- bind_rows(temp_asexual_ajuc, temp_sexual_ajuc)

ggplot(data = temp_df_ajuc, aes(x = reproductive_mode, y = temp_stability_scaled, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

overall_asexual_ajuc <- raster::extract(overall_stability_scaled, aste_locs_asexual) %>% 
  enframe(name = NULL, value = "overall_stability_scaled") %>% 
  mutate(reproductive_mode = "asexual")

overall_sexual_ajuc <- raster::extract(overall_stability_scaled, aste_locs_sexual) %>% 
  enframe(name = NULL, value = "overall_stability_scaled") %>% 
  mutate(reproductive_mode = "sexual")

overall_df_ajuc <- bind_rows(overall_asexual_ajuc, overall_sexual_ajuc)

ggplot(data = overall_df_ajuc, aes(x = reproductive_mode, y = overall_stability_scaled, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

evi_asexual_ajuc <- raster::extract(evi_stability_scaled, aste_locs_asexual) %>% 
  enframe(name = NULL, value = "evi_stability_scaled") %>% 
  mutate(reproductive_mode = "asexual")

evi_sexual_ajuc <- raster::extract(evi_stability_scaled, aste_locs_sexual) %>% 
  enframe(name = NULL, value = "evi_stability_scaled") %>% 
  mutate(reproductive_mode = "sexual")

evi_df_ajuc <- bind_rows(evi_asexual_ajuc, evi_sexual_ajuc)

ggplot(data = evi_df_ajuc, aes(x = reproductive_mode, y = evi_stability_scaled, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

Annual stability

evi_asexual_ajuc_annual <- raster::extract(evi_stability_scaled_annual, aste_locs_asexual) %>% 
  enframe(name = NULL, value = "evi_stability_scaled_annual") %>% 
  mutate(reproductive_mode = "asexual")

evi_sexual_ajuc_annual <- raster::extract(evi_stability_scaled_annual, aste_locs_sexual) %>% 
  enframe(name = NULL, value = "evi_stability_scaled_annual") %>% 
  mutate(reproductive_mode = "sexual")

evi_df_ajuc_annual <- bind_rows(evi_asexual_ajuc_annual, evi_sexual_ajuc_annual)

ggplot(data = evi_df_ajuc_annual, aes(x = reproductive_mode, y = evi_stability_scaled_annual, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

Clitarchus

summary_list_clita <- species_pca_fun(loc_clim, "clitarchus")
clita_plot <- plot_clim_pca(summary_list_clita$loc_clim, summary_list_clita$summary_pca, factor = "reproductive_mode")
clita_plot

#if selfcontained = TRUE, you can remove the folder that gets added alongside the plot. 
#htmlwidgets::saveWidget(clita_plot, paste0(getwd(), "/plots/repro_mode_plots/clitarchus_pca.html"), selfcontained = TRUE)

#filter localities for the focal genus
clita_loc <- loc %>% 
  filter(genus == "clitarchus")
  
#use sourced plot_locs_leaflet script to plot localities
clita_map <- plot_locs_leaflet(clita_loc, "reproductive_mode")

clita_map

#in case I want to save the map somewhere
#mapview::mapshot(clita_map, url = paste0(getwd(), "/plots/repro_mode_plots/clita_map.html"), file = paste0(getwd(), "/plots/repro_mode_plots/clita_map.pdf"))

summary_list_clita$summary_pca
loadings_clita <- summary_list_clita$summary_pca$rotation
knitr::kable(round(loadings_clita[,1:3],3)) #Table of loading scores for the first 3 PCs. 

Now I’m going to to environmental niche factor analysis between sexual and asexual populations within the species.

#get background env't for the species
choo_bg_env <- bg_env_crop(clita_loc, 
                           species = "hookeri",
                           environment = w, 
                           buffer = 0.5)

#enfa for the sexual species
choo_sexual_enfa <- enfa_calc_fun(locs = clita_loc, 
                                  species = "hookeri", 
                                  reproductive_mode = "sexual", 
                                  mask_raster = choo_bg_env)

#enfa for the asexual species
choo_asexual_enfa <- enfa_calc_fun(locs = clita_loc, 
                                   species = "hookeri", 
                                   reproductive_mode = "asexual", 
                                   mask_raster = choo_bg_env)


#plot the marginality scores
marginality_lollipop(sex_marg = choo_sexual_enfa$m, 
                    asex_marg = choo_asexual_enfa$m,
                    full_species_name = "Clitarchus hookeri")

A couple different ways to visualize the environmental variation. 1) Scatterplot visualizations of marginality vs axis 1 of the specialization with the labels removed (they make things indiscernable). Red = occupied e-space. Gray = Background e-space. 2) ENFA histogram visualizations of marginality and specialization axes. 3) PCA of total e-space with colors corresponding to sexual vs. asexual populations.

### 1) ENFA scatterplot
#access the relevant values for plotting
choo_asexual_df <- choo_asexual_enfa$li %>% 
  as_tibble() %>% 
  bind_cols(pr = choo_asexual_enfa$pr)


choo_sexual_df <- choo_sexual_enfa$li %>% 
  as_tibble() %>% 
  bind_cols(pr = choo_sexual_enfa$pr)


#asexual
enfa_hex_plot(choo_asexual_df, marg = Mar, spec = Spe1, repro_mode = "Asexual")


#sexual
enfa_hex_plot(choo_sexual_df, marg = Mar, spec = Spe1, repro_mode = "Sexual")


### 2) ENFA histogram
#asexual
hist(choo_asexual_enfa)
title(main = "Asexual", adj = 0.7, line = -12)

#sexual
hist(choo_sexual_enfa)
title(main = "Sexual", adj = 0.7, line = -12)

### 3) PCA of total e-space
total_climate_pca_plot(bg_env = choo_bg_env, locs = clita_loc, genus = "Clitarchus", species = "hookeri")

We’re also interested in seeing if asexual populations live in more unstable climatic areas relative to sexual populations. Looks like this pattern is similar to the general cross-species pattern.

clita_locs_asexual <- clita_loc %>% 
  mutate(lat_long = str_c(latitude, longitude, sep = "_")) %>% 
  filter(reproductive_mode == "asexual",
         !duplicated(lat_long)) %>% 
  dplyr::select(longitude, latitude)

clita_locs_sexual <- clita_loc %>% 
  mutate(lat_long = str_c(latitude, longitude, sep = "_")) %>% 
  filter(reproductive_mode == "sexual",
         !duplicated(lat_long)) %>% 
  dplyr::select(longitude, latitude)

precip_asexual_choo <- raster::extract(precip_stability_scaled, clita_locs_asexual) %>% 
  enframe(name = NULL, value = "precip_stability_scaled") %>% 
  mutate(reproductive_mode = "asexual")

precip_sexual_choo <- raster::extract(precip_stability_scaled, clita_locs_sexual) %>% 
  enframe(name = NULL, value = "precip_stability_scaled") %>% 
  mutate(reproductive_mode = "sexual")

precip_df_choo <- bind_rows(precip_asexual_choo, precip_sexual_choo)

ggplot(data = precip_df_choo, aes(x = reproductive_mode, y = precip_stability_scaled, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()


temp_asexual_choo <- raster::extract(temp_stability_scaled, clita_locs_asexual) %>% 
  enframe(name = NULL, value = "temp_stability_scaled") %>% 
  mutate(reproductive_mode = "asexual")

temp_sexual_choo <- raster::extract(temp_stability_scaled, clita_locs_sexual) %>% 
  enframe(name = NULL, value = "temp_stability_scaled") %>% 
  mutate(reproductive_mode = "sexual")

temp_df_choo <- bind_rows(temp_asexual_choo, temp_sexual_choo)

ggplot(data = temp_df_choo, aes(x = reproductive_mode, y = temp_stability_scaled, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

overall_asexual_choo <- raster::extract(overall_stability_scaled, clita_locs_asexual) %>% 
  enframe(name = NULL, value = "overall_stability_scaled") %>% 
  mutate(reproductive_mode = "asexual")

overall_sexual_choo <- raster::extract(overall_stability_scaled, clita_locs_sexual) %>% 
  enframe(name = NULL, value = "overall_stability_scaled") %>% 
  mutate(reproductive_mode = "sexual")

overall_df_choo <- bind_rows(overall_asexual_choo, overall_sexual_choo)

ggplot(data = overall_df_choo, aes(x = reproductive_mode, y = overall_stability_scaled, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

evi_asexual_choo <- raster::extract(evi_stability_scaled, clita_locs_asexual) %>% 
  enframe(name = NULL, value = "evi_stability_scaled") %>% 
  mutate(reproductive_mode = "asexual")

evi_sexual_choo <- raster::extract(evi_stability_scaled, clita_locs_sexual) %>% 
  enframe(name = NULL, value = "evi_stability_scaled") %>% 
  mutate(reproductive_mode = "sexual")

evi_df_choo <- bind_rows(evi_asexual_choo, evi_sexual_choo)

ggplot(data = evi_df_choo, aes(x = reproductive_mode, y = evi_stability_scaled, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

Annual stability

evi_asexual_choo_annual <- raster::extract(evi_stability_scaled_annual, clita_locs_asexual) %>% 
  enframe(name = NULL, value = "evi_stability_scaled_annual") %>% 
  mutate(reproductive_mode = "asexual")

evi_sexual_choo_annual <- raster::extract(evi_stability_scaled_annual, clita_locs_sexual) %>% 
  enframe(name = NULL, value = "evi_stability_scaled_annual") %>% 
  mutate(reproductive_mode = "sexual")

evi_df_choo_annual <- bind_rows(evi_asexual_choo_annual, evi_sexual_choo_annual)

ggplot(data = evi_df_choo_annual, aes(x = reproductive_mode, y = evi_stability_scaled_annual, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

Micrarchus

summary_list_micra <- species_pca_fun(loc_clim, "micrarchus")
micra_plot <- plot_clim_pca(summary_list_micra$loc_clim, summary_list_micra$summary_pca, factor = "reproductive_mode")
micra_plot

#if selfcontained = TRUE, you can remove the folder that gets added alongside the plot. 
#htmlwidgets::saveWidget(micra_plot, paste0(getwd(), "/plots/repro_mode_plots/micrarchus_pca.html"), selfcontained = TRUE)

#filter localities for the focal genus
micra_loc <- loc %>% 
  filter(genus == "micrarchus")
  
#use sourced plot_locs_leaflet script to plot localities
micra_map <- plot_locs_leaflet(micra_loc, "reproductive_mode")

micra_map

#in case I want to save the map somewhere
#mapview::mapshot(micra_map, url = paste0(getwd(), "/plots/repro_mode_plots/micra_map.html"), file = paste0(getwd(), "/plots/repro_mode_plots/micra_map.pdf"))

summary_list_micra$summary_pca
loadings_micra <- summary_list_micra$summary_pca$rotation
knitr::kable(round(loadings_micra[,1:3],3)) #Table of loading scores for the first 3 PCs. 

Niveaphasma

summary_list_nive <- species_pca_fun(loc_clim, "niveaphasma")
nive_plot <- plot_clim_pca(summary_list_nive$loc_clim, summary_list_nive$summary_pca, factor = "reproductive_mode")

nive_plot

#if selfcontained = TRUE, you can remove the folder that gets added alongside the plot. 
#htmlwidgets::saveWidget(nive_plot, paste0(getwd(), "/plots/repro_mode_plots/niveaphasma_pca.html"), selfcontained = TRUE)

#filter localities for the focal genus
nive_loc <- loc %>% 
  filter(genus == "niveaphasma")
  
#use sourced plot_locs_leaflet script to plot localities
nive_map <- plot_locs_leaflet(nive_loc, "reproductive_mode")

nive_map

#in case I want to save the map somewhere
#mapview::mapshot(nive_map, url = paste0(getwd(), "/plots/repro_mode_plots/nive_map.html"), file = paste0(getwd(), "/plots/repro_mode_plots/nive_map.pdf"))

summary_list_nive$summary_pca
loadings_nive <- summary_list_nive$summary_pca$rotation
knitr::kable(round(loadings_nive[,1:3],3)) #Table of loading scores for the first 3 PCs. 

Now I’m going to to environmental niche factor analysis between sexual and asexual populations within the species.

#get background env't for the species
nive_bg_env <- bg_env_crop(nive_loc, 
                           species = "annulata",
                           environment = w, 
                           buffer = 0.5)

#enfa for the sexual species
nive_sexual_enfa <- enfa_calc_fun(locs = nive_loc, 
                                  species = "annulata", 
                                  reproductive_mode = "sexual", 
                                  mask_raster = nive_bg_env)

#enfa for the asexual species
nive_asexual_enfa <- enfa_calc_fun(locs = nive_loc, 
                                   species = "annulata", 
                                   reproductive_mode = "asexual", 
                                   mask_raster = nive_bg_env)


#plot the marginality scores
marginality_lollipop(sex_marg = nive_sexual_enfa$m, 
                    asex_marg = nive_asexual_enfa$m,
                    full_species_name = "Niveaphasma annulata")

A couple different ways to visualize the environmental variation. 1) Scatterplot visualizations of marginality vs axis 1 of the specialization with the labels removed (they make things indiscernable). Red = occupied e-space. Gray = Background e-space. 2) ENFA histogram visualizations of marginality and specialization axes. 3) PCA of total e-space with colors corresponding to sexual vs. asexual populations.

### 1) ENFA scatterplot
#access the relevant values for plotting
nive_asexual_df <- nive_asexual_enfa$li %>% 
  as_tibble() %>% 
  bind_cols(pr = nive_asexual_enfa$pr)


nive_sexual_df <- nive_sexual_enfa$li %>% 
  as_tibble() %>% 
  bind_cols(pr = nive_sexual_enfa$pr)


#asexual
enfa_hex_plot(nive_asexual_df, marg = Mar, spec = Spe1, repro_mode = "Asexual")


#sexual
enfa_hex_plot(nive_sexual_df, marg = Mar, spec = Spe1, repro_mode = "Sexual")


### 2) ENFA histogram
#asexual
hist(nive_asexual_enfa)
title(main = "Asexual", adj = 0.7, line = -12)

#sexual
hist(nive_sexual_enfa)
title(main = "Sexual", adj = 0.7, line = -12)

### 3) PCA of total e-space
total_climate_pca_plot(bg_env = nive_bg_env, locs = nive_loc, genus = "Niveaphasma", species = "annulata")

We’re also interested in seeing if asexual populations live in more unstable climatic areas relative to sexual populations.

nive_locs_asexual <- nive_loc %>% 
  mutate(lat_long = str_c(latitude, longitude, sep = "_")) %>% 
  filter(reproductive_mode == "asexual",
         !duplicated(lat_long)) %>% 
  dplyr::select(longitude, latitude)

nive_locs_sexual <- nive_loc %>% 
  mutate(lat_long = str_c(latitude, longitude, sep = "_")) %>% 
  filter(reproductive_mode == "sexual",
         !duplicated(lat_long)) %>% 
  dplyr::select(longitude, latitude)

precip_asexual_nive <- raster::extract(precip_stability_scaled, nive_locs_asexual) %>% 
  enframe(name = NULL, value = "precip_stability_scaled") %>% 
  mutate(reproductive_mode = "asexual")

precip_sexual_nive <- raster::extract(precip_stability_scaled, nive_locs_sexual) %>% 
  enframe(name = NULL, value = "precip_stability_scaled") %>% 
  mutate(reproductive_mode = "sexual")

precip_df_nive <- bind_rows(precip_asexual_nive, precip_sexual_nive)

ggplot(data = precip_df_nive, aes(x = reproductive_mode, y = precip_stability_scaled, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()


temp_asexual_nive <- raster::extract(temp_stability_scaled, nive_locs_asexual) %>% 
  enframe(name = NULL, value = "temp_stability_scaled") %>% 
  mutate(reproductive_mode = "asexual")

temp_sexual_nive <- raster::extract(temp_stability_scaled, nive_locs_sexual) %>% 
  enframe(name = NULL, value = "temp_stability_scaled") %>% 
  mutate(reproductive_mode = "sexual")

temp_df_nive <- bind_rows(temp_asexual_nive, temp_sexual_nive)

ggplot(data = temp_df_nive, aes(x = reproductive_mode, y = temp_stability_scaled, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

overall_asexual_nive <- raster::extract(overall_stability_scaled, nive_locs_asexual) %>% 
  enframe(name = NULL, value = "overall_stability_scaled") %>% 
  mutate(reproductive_mode = "asexual")

overall_sexual_nive <- raster::extract(overall_stability_scaled, nive_locs_sexual) %>% 
  enframe(name = NULL, value = "overall_stability_scaled") %>% 
  mutate(reproductive_mode = "sexual")

overall_df_nive <- bind_rows(overall_asexual_nive, overall_sexual_nive)

ggplot(data = overall_df_nive, aes(x = reproductive_mode, y = overall_stability_scaled, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

evi_asexual_nive <- raster::extract(evi_stability_scaled, nive_locs_asexual) %>% 
  enframe(name = NULL, value = "evi_stability_scaled") %>% 
  mutate(reproductive_mode = "asexual")

evi_sexual_nive <- raster::extract(evi_stability_scaled, nive_locs_sexual) %>% 
  enframe(name = NULL, value = "evi_stability_scaled") %>% 
  mutate(reproductive_mode = "sexual")

evi_df_nive <- bind_rows(evi_asexual_nive, evi_sexual_nive)

ggplot(data = evi_df_nive, aes(x = reproductive_mode, y = evi_stability_scaled, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

Annual stability

evi_asexual_nive_annual <- raster::extract(evi_stability_scaled_annual, nive_locs_asexual) %>% 
  enframe(name = NULL, value = "evi_stability_scaled_annual") %>% 
  mutate(reproductive_mode = "asexual")

evi_sexual_nive_annual <- raster::extract(evi_stability_scaled_annual, nive_locs_sexual) %>% 
  enframe(name = NULL, value = "evi_stability_scaled_annual") %>% 
  mutate(reproductive_mode = "sexual")

evi_df_nive_annual <- bind_rows(evi_asexual_nive_annual, evi_sexual_nive_annual)

ggplot(data = evi_df_nive_annual, aes(x = reproductive_mode, y = evi_stability_scaled_annual, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

Spinotectarchus

summary_list_spin <- species_pca_fun(loc_clim, "spinotectarchus")
spin_plot <- plot_clim_pca(summary_list_spin$loc_clim, summary_list_spin$summary_pca, factor = "reproductive_mode")
spin_plot

#if selfcontained = TRUE, you can remove the folder that gets added alongside the plot. 
#htmlwidgets::saveWidget(spin_plot, paste0(getwd(), "/plots/repro_mode_plots/spinotectarchus_pca.html"), selfcontained = TRUE)

#filter localities for the focal genus
spin_loc <- loc %>% 
  filter(genus == "spinotectarchus")
  
#use sourced plot_locs_leaflet script to plot localities
spin_map <- plot_locs_leaflet(spin_loc, "reproductive_mode")

spin_map

#in case I want to save the map somewhere
#mapview::mapshot(spin_map, url = paste0(getwd(), "/plots/repro_mode_plots/spin_map.html"), file = paste0(getwd(), "/plots/repro_mode_plots/spin_map.pdf"))

summary_list_spin$summary_pca
loadings_spin <- summary_list_spin$summary_pca$rotation
knitr::kable(round(loadings_spin[,1:3],3)) #Table of loading scores for the first 3 PCs. 

Tectarchus

summary_list_tect <- species_pca_fun(loc_clim, "tectarchus")
tect_plot <- plot_clim_pca(summary_list_tect$loc_clim, summary_list_tect$summary_pca, factor = "reproductive_mode")
tect_plot

#if selfcontained = TRUE, you can remove the folder that gets added alongside the plot. 
#htmlwidgets::saveWidget(tect_plot, paste0(getwd(), "/plots/repro_mode_plots/tectarchus_pca.html"), selfcontained = TRUE)

#filter localities for the focal genus
tect_loc <- loc %>% 
  filter(genus == "tectarchus")
  
#use sourced plot_locs_leaflet script to plot localities
tect_map <- plot_locs_leaflet(tect_loc, "reproductive_mode")

tect_map

#in case I want to save the map somewhere
#mapview::mapshot(tect_map, url = paste0(getwd(), "/plots/repro_mode_plots/tect_map.html"), file = paste0(getwd(), "/plots/repro_mode_plots/tect_map.pdf"))

summary_list_tect$summary_pca
loadings_tect <- summary_list_tect$summary_pca$rotation
knitr::kable(round(loadings_tect[,1:3],3)) #Table of loading scores for the first 3 PCs. 

Now I’m going to to environmental niche factor analysis between sexual and asexual populations within the species.

This is for Tectarchus ovobessus.

#get background env't for the species
tect_ovo_bg_env <- bg_env_crop(tect_loc, 
                           species = "ovobessus",
                           environment = w, 
                           buffer = 0.5)

#enfa for the sexual species
tect_ovo_sexual_enfa <- enfa_calc_fun(locs = tect_loc, 
                                  species = "ovobessus", 
                                  reproductive_mode = "sexual", 
                                  mask_raster = tect_ovo_bg_env)

#enfa for the asexual species
tect_ovo_asexual_enfa <- enfa_calc_fun(locs = tect_loc, 
                                   species = "ovobessus", 
                                   reproductive_mode = "asexual", 
                                   mask_raster = tect_ovo_bg_env)


#plot the marginality scores
marginality_lollipop(sex_marg = tect_ovo_sexual_enfa$m, 
                    asex_marg = tect_ovo_asexual_enfa$m,
                    full_species_name = "Tectarchus ovobessus")

A couple different ways to visualize the environmental variation. 1) Scatterplot visualizations of marginality vs axis 1 of the specialization with the labels removed (they make things indiscernable). Red = occupied e-space. Gray = Background e-space. 2) ENFA histogram visualizations of marginality and specialization axes. 3) PCA of total e-space with colors corresponding to sexual vs. asexual populations.

### 1) ENFA scatterplot
#access the relevant values for plotting
tect_ovo_asexual_df <- tect_ovo_asexual_enfa$li %>% 
  as_tibble() %>% 
  bind_cols(pr = tect_ovo_asexual_enfa$pr)


tect_ovo_sexual_df <- tect_ovo_sexual_enfa$li %>% 
  as_tibble() %>% 
  bind_cols(pr = tect_ovo_sexual_enfa$pr)


#asexual
enfa_hex_plot(tect_ovo_asexual_df, marg = Mar, spec = Spe1, repro_mode = "Asexual")


#sexual
enfa_hex_plot(tect_ovo_sexual_df, marg = Mar, spec = Spe1, repro_mode = "Sexual")


### 2) ENFA histogram
#asexual
hist(tect_ovo_asexual_enfa)
title(main = "Asexual", adj = 0.7, line = -12)

#sexual
hist(tect_ovo_sexual_enfa)
title(main = "Sexual", adj = 0.7, line = -12)

### 3) PCA of total e-space
total_climate_pca_plot(bg_env = tect_ovo_bg_env, locs = tect_loc, genus = "Tectarchus", species = "ovobessus")

We’re also interested in seeing if asexual populations live in more unstable climatic areas relative to sexual populations.

tect_ovo_locs_asexual <- tect_loc %>% 
  mutate(lat_long = str_c(latitude, longitude, sep = "_")) %>% 
  filter(species == "ovobessus", 
         reproductive_mode == "asexual",
         !duplicated(lat_long)) %>% 
  dplyr::select(longitude, latitude)

tect_ovo_locs_sexual <- tect_loc %>% 
  mutate(lat_long = str_c(latitude, longitude, sep = "_")) %>% 
  filter(species == "ovobessus", 
         reproductive_mode == "sexual",
         !duplicated(lat_long)) %>% 
  dplyr::select(longitude, latitude)

precip_asexual_tect_ovo <- raster::extract(precip_stability_scaled, tect_ovo_locs_asexual) %>% 
  enframe(name = NULL, value = "precip_stability_scaled") %>% 
  mutate(reproductive_mode = "asexual")

precip_sexual_tect_ovo <- raster::extract(precip_stability_scaled, tect_ovo_locs_sexual) %>% 
  enframe(name = NULL, value = "precip_stability_scaled") %>% 
  mutate(reproductive_mode = "sexual")

precip_df_tect_ovo <- bind_rows(precip_asexual_tect_ovo, precip_sexual_tect_ovo)

ggplot(data = precip_df_tect_ovo, aes(x = reproductive_mode, y = precip_stability_scaled, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()


temp_asexual_tect_ovo <- raster::extract(temp_stability_scaled, tect_ovo_locs_asexual) %>% 
  enframe(name = NULL, value = "temp_stability_scaled") %>% 
  mutate(reproductive_mode = "asexual")

temp_sexual_tect_ovo <- raster::extract(temp_stability_scaled, tect_ovo_locs_sexual) %>% 
  enframe(name = NULL, value = "temp_stability_scaled") %>% 
  mutate(reproductive_mode = "sexual")

temp_df_tect_ovo <- bind_rows(temp_asexual_tect_ovo, temp_sexual_tect_ovo)

ggplot(data = temp_df_tect_ovo, aes(x = reproductive_mode, y = temp_stability_scaled, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

overall_asexual_tect_ovo <- raster::extract(overall_stability_scaled, tect_ovo_locs_asexual) %>% 
  enframe(name = NULL, value = "overall_stability_scaled") %>% 
  mutate(reproductive_mode = "asexual")

overall_sexual_tect_ovo <- raster::extract(overall_stability_scaled, tect_ovo_locs_sexual) %>% 
  enframe(name = NULL, value = "overall_stability_scaled") %>% 
  mutate(reproductive_mode = "sexual")

overall_df_tect_ovo <- bind_rows(overall_asexual_tect_ovo, overall_sexual_tect_ovo)

ggplot(data = overall_df_tect_ovo, aes(x = reproductive_mode, y = overall_stability_scaled, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

evi_asexual_tect_ovo <- raster::extract(evi_stability_scaled, tect_ovo_locs_asexual) %>% 
  enframe(name = NULL, value = "evi_stability_scaled") %>% 
  mutate(reproductive_mode = "asexual")

evi_sexual_tect_ovo <- raster::extract(evi_stability_scaled, tect_ovo_locs_sexual) %>% 
  enframe(name = NULL, value = "evi_stability_scaled") %>% 
  mutate(reproductive_mode = "sexual")

evi_df_tect_ovo <- bind_rows(evi_asexual_tect_ovo, evi_sexual_tect_ovo)

ggplot(data = evi_df_tect_ovo, aes(x = reproductive_mode, y = evi_stability_scaled, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

Annual stability

evi_asexual_tect_ovo_annual <- raster::extract(evi_stability_scaled_annual, tect_ovo_locs_asexual) %>% 
  enframe(name = NULL, value = "evi_stability_scaled_annual") %>% 
  mutate(reproductive_mode = "asexual")

evi_sexual_tect_ovo_annual <- raster::extract(evi_stability_scaled_annual, tect_ovo_locs_sexual) %>% 
  enframe(name = NULL, value = "evi_stability_scaled_annual") %>% 
  mutate(reproductive_mode = "sexual")

evi_df_tect_ovo_annual <- bind_rows(evi_asexual_tect_ovo_annual, evi_sexual_tect_ovo_annual)

ggplot(data = evi_df_tect_ovo_annual, aes(x = reproductive_mode, y = evi_stability_scaled_annual, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

This is an enfa for Tectarchus huttoni.

###Only need to get bg env't if you're not running the previous chunk
#get background env't for the species
tect_hutt_bg_env <- bg_env_crop(tect_loc, 
                           species = "huttoni",
                           environment = w, 
                           buffer = 0.5)

#enfa for the sexual species
tect_hutt_sexual_enfa <- enfa_calc_fun(locs = tect_loc, 
                                  species = "huttoni", 
                                  reproductive_mode = "sexual", 
                                  mask_raster = tect_hutt_bg_env)

#enfa for the asexual species
tect_hutt_asexual_enfa <- enfa_calc_fun(locs = tect_loc, 
                                   species = "huttoni", 
                                   reproductive_mode = "asexual", 
                                   mask_raster = tect_hutt_bg_env)


#plot the marginality scores
marginality_lollipop(sex_marg = tect_hutt_sexual_enfa$m, 
                    asex_marg = tect_hutt_asexual_enfa$m,
                    full_species_name = "Tectarchus huttoni")

A couple different ways to visualize the environmental variation. 1) Scatterplot visualizations of marginality vs axis 1 of the specialization with the labels removed (they make things indiscernable). Red = occupied e-space. Gray = Background e-space. 2) ENFA histogram visualizations of marginality and specialization axes. 3) PCA of total e-space with colors corresponding to sexual vs. asexual populations.

### 1) ENFA scatterplot
#access the relevant values for plotting
tect_hutt_asexual_df <- tect_hutt_asexual_enfa$li %>% 
  as_tibble() %>% 
  bind_cols(pr = tect_hutt_asexual_enfa$pr)


tect_hutt_sexual_df <- tect_hutt_sexual_enfa$li %>% 
  as_tibble() %>% 
  bind_cols(pr = tect_hutt_sexual_enfa$pr)


#asexual
enfa_hex_plot(tect_hutt_asexual_df, marg = Mar, spec = Spe1, repro_mode = "Asexual")


#sexual
enfa_hex_plot(tect_hutt_sexual_df, marg = Mar, spec = Spe1, repro_mode = "Sexual")


### 2) ENFA histogram
#asexual
hist(tect_hutt_asexual_enfa)
title(main = "Asexual", adj = 0.7, line = -12)

#sexual
hist(tect_hutt_sexual_enfa)
title(main = "Sexual", adj = 0.7, line = -12)

### 3) PCA of total e-space
total_climate_pca_plot(bg_env = tect_hutt_bg_env, locs = tect_loc, genus = "Tectarchus", species = "huttoni")

tect_hutt_locs_asexual <- tect_loc %>% 
  mutate(lat_long = str_c(latitude, longitude, sep = "_")) %>% 
  filter(species == "huttoni", 
         reproductive_mode == "asexual",
         !duplicated(lat_long)) %>% 
  dplyr::select(longitude, latitude)

tect_hutt_locs_sexual <- tect_loc %>% 
  mutate(lat_long = str_c(latitude, longitude, sep = "_")) %>% 
  filter(species == "huttoni", 
         reproductive_mode == "sexual",
         !duplicated(lat_long)) %>% 
  dplyr::select(longitude, latitude)

precip_asexual_tect_hutt <- raster::extract(precip_stability_scaled, tect_hutt_locs_asexual) %>% 
  enframe(name = NULL, value = "precip_stability_scaled") %>% 
  mutate(reproductive_mode = "asexual")

precip_sexual_tect_hutt <- raster::extract(precip_stability_scaled, tect_hutt_locs_sexual) %>% 
  enframe(name = NULL, value = "precip_stability_scaled") %>% 
  mutate(reproductive_mode = "sexual")

precip_df_tect_hutt <- bind_rows(precip_asexual_tect_hutt, precip_sexual_tect_hutt)

ggplot(data = precip_df_tect_hutt, aes(x = reproductive_mode, y = precip_stability_scaled, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()


temp_asexual_tect_hutt <- raster::extract(temp_stability_scaled, tect_hutt_locs_asexual) %>% 
  enframe(name = NULL, value = "temp_stability_scaled") %>% 
  mutate(reproductive_mode = "asexual")

temp_sexual_tect_hutt <- raster::extract(temp_stability_scaled, tect_hutt_locs_sexual) %>% 
  enframe(name = NULL, value = "temp_stability_scaled") %>% 
  mutate(reproductive_mode = "sexual")

temp_df_tect_hutt <- bind_rows(temp_asexual_tect_hutt, temp_sexual_tect_hutt)

ggplot(data = temp_df_tect_hutt, aes(x = reproductive_mode, y = temp_stability_scaled, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

overall_asexual_tect_hutt <- raster::extract(overall_stability_scaled, tect_hutt_locs_asexual) %>% 
  enframe(name = NULL, value = "overall_stability_scaled") %>% 
  mutate(reproductive_mode = "asexual")

overall_sexual_tect_hutt <- raster::extract(overall_stability_scaled, tect_hutt_locs_sexual) %>% 
  enframe(name = NULL, value = "overall_stability_scaled") %>% 
  mutate(reproductive_mode = "sexual")

overall_df_tect_hutt <- bind_rows(overall_asexual_tect_hutt, overall_sexual_tect_hutt)

ggplot(data = overall_df_tect_hutt, aes(x = reproductive_mode, y = overall_stability_scaled, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

evi_asexual_tect_hutt <- raster::extract(evi_stability_scaled, tect_hutt_locs_asexual) %>% 
  enframe(name = NULL, value = "evi_stability_scaled") %>% 
  mutate(reproductive_mode = "asexual")

evi_sexual_tect_hutt <- raster::extract(evi_stability_scaled, tect_hutt_locs_sexual) %>% 
  enframe(name = NULL, value = "evi_stability_scaled") %>% 
  mutate(reproductive_mode = "sexual")

evi_df_tect_hutt <- bind_rows(evi_asexual_tect_hutt, evi_sexual_tect_hutt)

ggplot(data = evi_df_tect_hutt, aes(x = reproductive_mode, y = evi_stability_scaled, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

Annual stability

evi_asexual_tect_hutt_annual <- raster::extract(evi_stability_scaled_annual, tect_hutt_locs_asexual) %>% 
  enframe(name = NULL, value = "evi_stability_scaled_annual") %>% 
  mutate(reproductive_mode = "asexual")

evi_sexual_tect_hutt_annual <- raster::extract(evi_stability_scaled_annual, tect_hutt_locs_sexual) %>% 
  enframe(name = NULL, value = "evi_stability_scaled_annual") %>% 
  mutate(reproductive_mode = "sexual")

evi_df_tect_hutt_annual <- bind_rows(evi_asexual_tect_hutt_annual, evi_sexual_tect_hutt_annual)

ggplot(data = evi_df_tect_hutt_annual, aes(x = reproductive_mode, y = evi_stability_scaled_annual, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

Tepakiphasma

Nothing. Only one locality.

Convenience scripts

These scripts aren’t crucial for reproducing this analysis, but were helpful for various reasons. Some of these have hard-coded paths and such, so no guarantees for use right out of the box.

This was a script I used to take full chelsa files, crop them to New Zealand extent, and write them to a file with my personal computer. I don’t have much memory, so unzipping to a temporary directory, then deleting the directory to free up space for the large files worked.

## get chelsa data
chelsa_folder <- "/Users/connorfrench/Dropbox/Old_Mac/climate-data/chelsa_30s_bio"
zip_files <- list.files(chelsa_folder, full.names = TRUE)

# using the Unarchiver commandline tools for Mac to unzip the 7zip chelsa layers. Regular unzip() does not work with 7z zipped files
for (file in zip_files) {
  # set temp directory
  tempd <- tempdir()
  system(paste("unar", file, "-o", tempd))
  r <- raster(list.files(tempd, pattern = "*.tif", full.names = TRUE)) %>%
    crop(extent(166, 179, -48, -34))
  writeRaster(r, filename = paste0("~/Desktop/", list.files(tempd, pattern = "*.tif")), format = "GTiff")
  unlink(tempd, recursive = TRUE)
}